 /* test the validity of the series following "Number Counts and Non-Gaussianity"
  * arXiv: 1311.7361
  */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "series.h"

int main(int argc, char *argv[]) {
  // compute the errors for different values of M3 and nuc
  double M3i, NUCi, nucxmax, nucx, eror, M3s;
  double temp, errtol;
  int N=atoi(argv[1]);
  FILE *fp;
  char fname[200];
  int feeder=0;
<<<<<<< .mine
  M3s=0.01;
=======
  M3s=0.1;
>>>>>>> .r23
  if (argv[2]) {
    feeder = atoi(argv[2]);
<<<<<<< .mine
    if (feeder) M3s=0.001;
=======
    if (feeder) M3s=0.01;
>>>>>>> .r23
  }

  for (M3i=M3s; M3i<M3s*12; M3i=M3i+M3s) {
    sprintf(fname, "/home/saroj/Research/TwoField/Simulation/Theory/sdata/series%d_%g_%d.txt",N, M3i, feeder);
    fp = fopen(fname, "w");
    for (NUCi=0.001; NUCi<10; NUCi=NUCi+0.001) {
      // look in the range of Nuci<nu<numax and find the max error percentange
      eror=pnextf(N, NUCi, M3i, feeder)/pnf(N, NUCi, M3i, feeder);
	//if (temp==0) printf("N: %d, nucx: %g, M3: %g, f: %d\n", N, nucx, M3i, feeder)
      
      fprintf(fp, "%g %g\n", NUCi, eror);
      //printf("M3i: %g, NUCi %g, error: %g\n", M3i, NUCi, error);
    }
    fclose(fp);
  } 

 // work in increments of M3s/2 and check for 20% error
/*  
  M3i=M3s;
  errtol=0.1;
  sprintf(fname, "/home/saroj/Research/feedersim/Theory/sdata/m3max%d_%d_%g.txt", N, feeder,errtol);
  fp = fopen(fname, "w");
  printf("oen\n");
<<<<<<< .mine
  for (NUCi=-1; NUCi<10; NUCi=NUCi+0.05) {
=======
  for (NUCi=0.01; NUCi<10; NUCi=NUCi+0.1) {
>>>>>>> .r23
    nucxmax=numaxf(NUCi, feeder);
    eror=0.;    
    while (fabs(eror)<errtol) {
      for (nucx=NUCi; nucx<nucxmax; nucx=nucx+0.01) {   
        temp=pnextf(N, nucx, M3i, feeder)/pnf(N, nucx, M3i, feeder);
        if (fabs(temp)>fabs(eror)) {
          eror=temp;
        }
      }
      M3i=M3i+M3s;
	  }
	  fprintf(fp, "%g %g\n", NUCi, M3i);
	  M3i=M3s;
	}
	fclose(fp);

  

  // check if the procedure works:
  printf("Hen(9,1): %g, direct: %g\n", Hen(8, .1), He8(.1));
  //printf("%g\n",pif(16, 4, .01, 1));
 // printf("pnextf(6, 1, .1): %g, direct: %g \n", pnextf(16, 4, .01, 1), p3f(6, .1));
  pif1(N, 2, .1);
  printf("next\n");
  pif1(N+1, 2, .1);
  //check
  printf("p1f1(4, 2, .1): %g, p4f(2, .1): %g\n", pif1(7, 2, .1), p7f(2,.1));
*/
}



double pnf(int n, double nu, double M3, int feeder) {
  double val=1.;
  int i;
  for (i=1; i<=n; i++) {
    if (feeder) {
      val+=pif1(i, nu, M3);
    } else {
      val+=pif(i, nu, M3);
    }
  }
  return val;
  
}

double pnextf(int n, double nu, double M3, int feeder) {
//  if (n==7) {return p8f(nu, M3);}
  if (feeder) {
    return pif1(n+1, nu, M3);
  } else {
    return pif(n+1, nu, M3);
  }
}

double numaxf(double nuc, int feeder) {
  if (feeder) {
    return 2.1 * pow(nuc, 0.7); // from 1311.7361 section 2.2
  } else {
    return 2.2 * pow(nuc, 0.7);
  }
}


double pif1(int i, double nu, double M3) {
  int nsol=0; // # of solns for this term
  i=i+2;
  int k[i+1]; // temporarily store the k's i.e the solns
  int j;
  int mold=1;
  int m, sumcur, r;
  double sum=0., prod=1.;
  for (j=1; j<=i; j++) {k[j]=0;} k[1]=i; // this is one known solutions
  //sum+=Hen(i+2*i, nu)*pow(M3/6, i)/tgamma(i+1.); // fix this
  if (i==3) {sum+=Hen(3, nu)*M3/6.;}
  while (mold!=i) {
    m = 1;
    sumcur= i;
  
    do {
      sumcur = sumcur - k[m]*m+m+1;
      //printf("m: %d , sumcur: %d\n", m, sumcur);
      k[m]=0;
      k[m+1]=k[m+1]+1;
      m=m+1;
    } while (sumcur > i && m <= mold);
    if (m>mold) {mold = m;}
    k[1]=i-sumcur;

    if (k[1]==0 && k[2]==0) {
    nsol=nsol+1; //printf("\n"); for (j=1; j<=i;j++) {printf("%d ", k[j]);};
    // sum over these k's to get the ith term p_i of the (h) pdf
    r=0; prod=1.0;
    for (j=1; j<=i-2; j++) {
      r += k[j];
      prod*=pow(Mnf(j+2, M3)/tgamma(j+3.),k[j+2])/tgamma(k[j+2]+1.);
    }
    sum+=prod*Hen(i, nu);
    //for (j=1; j<=i; j++) {printf("%d ", k[j]);} printf("\n %d \n",i );
    }
  }
  //printf("nsol: %d\n", nsol);
  return sum;
}  

double pif(int i, double nu, double M3) {
  // sum_{{k_m}} He_{i+2r}(nu) ||_{m=1}^i (M_{m+2}/(m+2)!)^{k_m}
  /* the algorithm to find the set {k_m} is given in "Expansions for nearly Gaussian expansions". {k_m} consists of all non-negative integer solutions of the Diophantine equation: k1 + 2 k2 + ... + n kn = n
  For example: for s=1, k1 = 1 so {1}
	       for s=2, k1 + 2 k2 = 2 so {0, 1}, {2, 0}
	       for s=3, k1 + 2 k2 + 3 k3 = 3 so {0, 0, 1}, {1, 1, 0}, {3,0 ,0}
	       and so on
	       
  */
  
  int nsol=1; // # of solns for this term
  int k[i+1]; // temporarily store the k's i.e the solns
  int j;
  int mold=1;
  int m, sumcur, r;
  double sum=0., prod=1.;
  for (j=1; j<=i; j++) {k[j]=0;} k[1]=i; // this is one known solutions
  sum+=Hen(i+2*i, nu)*pow(M3/6, i)/tgamma(i+1.);
  while (mold!=i) {
    m = 1;
    sumcur= i;
  
    do {
      sumcur = sumcur - k[m]*m+m+1;
      //printf("m: %d , sumcur: %d\n", m, sumcur);
      k[m]=0;
      k[m+1]=k[m+1]+1;
      m=m+1;
    } while (sumcur > i && m <= mold);
    if (m>mold) {mold = m;}
    k[1]=i-sumcur;
    nsol=nsol+1; //printf("\n"); for (j=1; j<=i;j++) {printf("%d ", k[j]);};
    // sum over these k's to get the ith term p_i of the (h) pdf
    r=0; prod=1.0;
    for (j=1; j<=i; j++) {
      r += k[j];
      prod*=pow(Mnh(j+2, M3)/tgamma(j+3.), k[j])/tgamma(k[j]+1.);
    }
    sum+=prod*Hen(i+2*r, nu);
  }
  //printf("nsold: %d\n", nsol);
  return sum;
}



double An(int n) {return tgamma(n+1.)*pow(2, n-3);}
double Bn(int n) {return tgamma(n)*pow(2, n-1);}
double Mnf(int n, double M3) {return Bn(n)*pow(M3/8., n/3);}
double Mnh(int n, double M3) {return An(n)*pow(M3/6., n-2);}

// the Hermite polynomials
double Hen(int n, double x) {
  // return the probabilits' Hermite polynomials directly for smal n and by computing the coefficients  from the recursion relation for n>9
  if (n==1) return He1(x); if (n==2) return He2(x); if (n==3) return He3(x);
  //if (n==4) return He4(x); if (n==5) return He5(x); if (n==6) return He6(x);
  //if (n==7) return He7(x); if (n==8) return He8(x); if (n==9) return He9(x);
  // if not compute the coefficients
  int a[n*3][n*3+1]; // can compute upto n
  int i,k;
  for (i=0; i<n*3; i++) for (k=0; k<n*3+1; k++) a[i][k]=0;
  
  a[1][0]=0; a[1][1]=1; a[0][0]=1;
  //printf("%d \n", a[0][0]);
  for (i=2; i<n+1; i++) {
    a[i][0]=-(i-1)*a[i-2][0];
    //printf("%d ", a[i][0]);
    for (k=1; k<=i; k++) {
      a[i][k]=a[i-1][k-1]-(i-1)*a[i-2][k];
      //printf("%d ", a[i][k]);
    }
    //printf("\n");
  }
  double sum = 0.;
  for (k=0; k<=n; k++) {
    sum+=a[n][k]*pow(x,k);
  }
  return sum;
}

// exact formulae for checking the algorithms used above for all orders.
double He1(double x) {return x;}
double He2(double x) {return x*x-1;}
double He3(double x) {return pow(x,3)-3*x;}
double He4(double x) {return pow(x,4)-6*x*x+3;}
double He5(double x) {return pow(x,5)-10*pow(x,3)+15*x;}
double He6(double x) {return pow(x,6)-15*pow(x,4)+45*x*x-15;}
double He7(double x) {return pow(x,7)-21*pow(x,5)+105*pow(x,3)-105*x;}
double He8(double x) {return pow(x,8)-28*pow(x,6)+210*pow(x,4)-420*x*x+105;}
double He9(double x) {return pow(x,9)-36*pow(x,7)+378*pow(x,5)-1260*pow(x,3)+945*x;}



double p1f(double nu, double M3) {
  return He3(nu)*M3/6.;
}
double p2f(double nu, double M3) {
  return He4(nu)*Mnf(4, M3)/24;
}
double p3f(double nu, double M3) {
  return He5(nu)*Mnf(5, M3)/120;
}
double p4f(double nu, double M3) {
  return He6(nu)*(Mnf(6, M3)/720+0.5*pow(M3/6,2));
}
double p5f(double nu, double M3) {
  return He7(nu)*(M3*Mnf(4, M3)/(6*24)+Mnf(7, M3)/(7*720));
}
double p6f(double nu, double M3) {
  return He8(nu)*(Mnf(8, M3)/(8*7*720)+0.5*pow(Mnf(4, M3)/24,2)+M3*Mnf(5, M3)/(6*120));
}
double p7f(double nu, double M3) {
  return He9(nu)*(Mnf(9, M3)/(9*8*7*720)+pow(M3/6,3)/6+(M3*Mnf(6, M3))/(6*720)+(Mnf(4, M3)*Mnf(5, M3))/(120*24));
}
